Integration of multimodal cues does not alter mean but reduces variance in avian responses to predators: a systematic review and meta-analysis

1 Setting up

Code

#install.packages("pacman")
#pacman::p_load(devtools, tidyverse, metafor, patchwork, R.rsp, emmeans)

#devtools::install_github("daniel1noble/orchaRd", force = TRUE)

library(tidyverse)
library(here)
library(lme4)
library(orchaRd)
#library(gptstudio)
library(metafor)
library(patchwork)
library(alluvial)
library(ggalluvial)
library(easyalluvial)
library(ape)
library(clubSandwich)
library(emmeans)
library(MuMIn)
# making metafor talk to MuMIn
eval(metafor:::.MuMIn)
# install.packages("pak")
#pak::pak("MichelNivard/gptstudio")

2 Getting data loaded

Code

#dat_full <- read.csv(here("data/dat_04_04_2023.csv"))
dat_full <- read.csv(here("data/dat_28_06_2023.csv"))

# add phylogenetic tree - only topologies
# TODO? - we could get better tree from birdtree.org
# we can do 50 different trees as in 
# https://academic.oup.com/sysbio/article/68/4/632/5267840
tree_top <- read.tree(here("R/birds_MA.tre"))

# tree with branch lengths
tree <- compute.brlen(tree_top)
#plot(tree)
# turning it into a correlation matrix
cor_tree <- vcv(tree,corr=T)

3 Preparing data set

Code

# function for calculating variance
Vd_func <- function(d, n1, n2, design, r = 0.5){
  # independent design
  if(design == "among"){
    var <- (n1 + n2) / (n1*n2) + d^2 / (2 * (n1 + n2 - 2)) # variance
  } else { # dependent design
    var <- 2*(1-r) / n1 + d^2 / (2*(n1 - 1)) # variance
  }
  var # return variance
}

# getting Hedges' g - get small size bias corrected effect size
dat_full$SMD <- dat_full$d / (1 - 3/(4 * (dat_full$NTreat + dat_full$Ncontrol) - 9))

# flipping d 
dat_full$SMD <- dat_full$d*dat_full$Direction*dat_full$PredictedDirection


# calucating Vd
dat_full$Vd <- with(dat_full, pmap_dbl(list(SMD, NTreat, Ncontrol, Design), Vd_func))


# extra useful function
# function for getting mean and sd from median, quartiles and sample size
# get_mean_sd <- function(median, q1, q3, n){
#   sd <- (q3 - q1) / (2 * (qnorm((0.75 * n - 0.125) / (n + 0.25)))) # sd
#   mean <- (median + q1 + q3)/3 # mean
#   c(mean, sd)
# }


# observation id
dat_full$Obs_ID <- 1:nrow(dat_full)
dat_full$Phylo <- gsub(" ", "_", dat_full$FocalSpL)

# filtering very large variance and also very small sample size
dat_int <- dat_full %>% filter(Vd < 10 & Ncontrol > 2 & NTreat > 2)

#dim(dat_full)
#dim(dat_int)


# sorting out modality stuff
# creat - 1,2,3 modality - also easier classification A, O, V (AOV = L) 

dat_int %>% mutate(Treat_mod = case_when(Treatment == "A" ~ "A",
                                          Treatment == "AV" ~ "AV",
                                          Treatment == "AVG" ~ "AV",
                                          Treatment == "AVM" ~ "AV",
                                          Treatment == "L" ~ "AVO",
                                          Treatment == "O" ~ "O",
                                          Treatment == "OV" ~ "OV",
                                          Treatment == "V" ~ "V",
                                          Treatment == "VG" ~ "V",
                                          Treatment == "VM" ~ "V",
                                          Treatment == "VP" ~ "V"),
                    # into how many
                    Treat_No = case_when(Treatment == "A" ~ 1,
                                         Treatment == "AV" ~ 2,
                                         Treatment == "AVG" ~ 2,
                                         Treatment == "AVM" ~ 2,
                                         Treatment == "L" ~ 3,
                                         Treatment == "O" ~ 1,
                                         Treatment == "OV" ~ 2,
                                         Treatment == "V" ~ 1,
                                         Treatment == "VG" ~ 1,
                                         Treatment == "VM" ~ 1,
                                         Treatment == "VP" ~ 1),
                    # des it have some add-ons
                    Add_on = case_when(Treatment == "A" ~ "No",
                                         Treatment == "AV" ~ "No",
                                         Treatment == "AVG" ~ "Yes",
                                         Treatment == "AVM" ~ "Yes",
                                         Treatment == "L" ~ "No",
                                         Treatment == "O" ~ "No",
                                         Treatment == "OV" ~ "No",
                                         Treatment == "V" ~ "No",
                                         Treatment == "VG" ~ "Yes",
                                         Treatment == "VM" ~ "Yes",
                                         Treatment == "VP" ~ "Yes"),

                      ) -> dat

# creating data just for A, V, and AV 
dat_short <- dat %>% filter(Treat_mod == "A" | Treat_mod == "V" | Treat_mod == "AV")

# for add-on, we only need V and AV
dat_short_add <- dat %>% filter(Treat_mod == "AV" | Treat_mod == "V")


dat <- dat %>%
  mutate_if(is.character, as.factor)

#summary(dat)

4 Exploratory visualization

Code

# some data exploration
# dat %>% group_by(Treat_mod) %>% summarise(n = n())
# using a table to see overalps between Treat_mod and Type
tab <- table(dat_short$Treat_mod, dat_short$Type)

# visualise this table using alluvial plot?
# https://cran.r-project.org/web/packages/alluvial/vignettes/alluvial.html

# Treatment vs Type
dat_short %>% group_by(Treat_mod, Type) %>%
  summarise(n = n()) -> tab1
#alluvial(tab1[,1:2], freq = tab1$n)

# using ggaruvial
ggplot(tab1,
       aes(y = n,
           axis1 = Treat_mod,
           axis2 = Type)) +
  geom_alluvium(aes(fill = Treat_mod)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 6, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Treatment modality and trait type")

Code

# Type vs Sex

dat %>% group_by(Sex, Type) %>%
  summarise(n = n()) -> tab2
#alluvial(tab2[,1:2], freq = tab2$n)

# using ggaruvial
ggplot(tab2,
       aes(y = n,
           axis1 = Type,
           axis2 = Sex)) +
  geom_alluvium(aes(fill = Type)) +
  geom_stratum(alpha = 0.5) +
  geom_text(stat = "stratum", size = 6, aes(label = after_stat(stratum))) +
  theme(legend.position = "none") +
  theme(legend.position = "none",
        axis.text.x = element_blank()) + # remove x-axis labels
  ylab("Frequency") + 
  xlab("Trait type and sex")

Code

# TOOD - use easyalluvial to make one fig
# other ones - easyalluvial
#https://www.r-bloggers.com/2018/10/data-exploration-with-alluvial-plots-an-introduction-to-easyalluvial/
# using easyalluvial and alluvial_wide
#factor_cols <- dat %>% select_if(is.factor) %>% names()

#alluvial_wide(dat_short, , max_variables = 5
#                , fill_by = 'first_variable' ) %>%
#  add_marginal_histograms(dat_short)

# exploratory analysis
# check each columns for missing values and other stuff
# dat %>% map_df(~sum(is.na(.)))

5 Meta-analysis

5.1 All random effects

Code

##| warning: false
# VCV matrix

VCV <- vcalc(vi = dat$Vd,
             cluster = dat$SubjectID,
             rho = 0.5)

mod0 <- rma.mv(yi = SMD,
       V = VCV, 
       random = list(~1 | Phylo,
                     ~1 | FocalSpL,
                     ~1 | RecNo,
                     ~1 | SubjectID, # incoprated as VCV
                     ~1 | Obs_ID),
       R = list(Phylo = cor_tree),
       test = "t",
       method = "REML", 
       sparse = TRUE,
       data = dat)

summary(mod0)

Multivariate Meta-Analysis Model (k = 629; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-881.6503  1763.3005  1775.3005  1801.9558  1775.4358   

Variance Components:

            estim    sqrt  nlvls  fixed     factor    R 
sigma^2.1  0.0000  0.0000     87     no      Phylo  yes 
sigma^2.2  0.0201  0.1418     87     no   FocalSpL   no 
sigma^2.3  0.1495  0.3866    114     no      RecNo   no 
sigma^2.4  0.0000  0.0000    145     no  SubjectID   no 
sigma^2.5  0.7132  0.8445    629     no     Obs_ID   no 

Model Results:

estimate      se    tval   df    pval   ci.lb   ci.ub      
  0.4312  0.0690  6.2518  628  <.0001  0.2958  0.5667  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

# TODO - think about whether we add this or not
robust(mod0, cluster = dat$SubjectID)

Multivariate Meta-Analysis Model (k = 629; method: REML)

Variance Components:

            estim    sqrt  nlvls  fixed     factor    R 
sigma^2.1  0.0000  0.0000     87     no      Phylo  yes 
sigma^2.2  0.0201  0.1418     87     no   FocalSpL   no 
sigma^2.3  0.1495  0.3866    114     no      RecNo   no 
sigma^2.4  0.0000  0.0000    145     no  SubjectID   no 
sigma^2.5  0.7132  0.8445    629     no     Obs_ID   no 

Number of estimates:   629
Number of clusters:    145
Estimates per cluster: 1-29 (mean: 4.34, median: 2)

Model Results:

estimate      se¹    tval¹   df¹    pval¹   ci.lb¹   ci.ub¹      
  0.4312  0.0572   7.5421   144   <.0001   0.3182   0.5442   *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

1) results based on cluster-robust inference (var-cov estimator: CR1,
   approx t-test and confidence interval, df: residual method)

Code

round(i2_ml(mod0), 2)

    I2_Total     I2_Phylo  I2_FocalSpL     I2_RecNo I2_SubjectID    I2_Obs_ID 
       93.17         0.00         2.12        15.78         0.00        75.27 

Code

orchard_plot(mod0,
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

5.2 Reduced model

Code

# reduced model

mod0r <- rma.mv(yi = SMD,
       V = VCV, 
       random = list(#~1 | Phylo,
                     ~1 | FocalSpL,
                     ~1 | RecNo,
                     #~1 | SubjectID, # incoprated as VCV
                     ~1 | Obs_ID),
       #R = list(Phylo = cor_tree),
       test = "t",
       method = "REML", 
       sparse = TRUE,
       data = dat)

summary(mod0r)

Multivariate Meta-Analysis Model (k = 629; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-881.6503  1763.3005  1771.3005  1789.0707  1771.3647   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0201  0.1418     87     no  FocalSpL 
sigma^2.2  0.1495  0.3866    114     no     RecNo 
sigma^2.3  0.7132  0.8445    629     no    Obs_ID 

Model Results:

estimate      se    tval   df    pval   ci.lb   ci.ub      
  0.4312  0.0690  6.2518  628  <.0001  0.2958  0.5667  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(i2_ml(mod0r), 2)

   I2_Total I2_FocalSpL    I2_RecNo   I2_Obs_ID 
      93.17        2.12       15.78       75.27 

6 Meta-regression: uni-moderator (mostly)

6.1 Treatmeant with all data

Code

## Treatment - A, V, AV etc 
#dat$Phylo <- as.character(dat$Phylo)
#match(dat$Phylo, colnames(cor_tree))
#match(colnames(cor_tree), dat$Phylo)

mod1 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1 | FocalSpL,
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mod = ~ Treat_mod - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)

summary(mod1)

Multivariate Meta-Analysis Model (k = 629; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-874.5965  1749.1931  1767.1931  1807.1040  1767.4867   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0135  0.1161     87     no  FocalSpL 
sigma^2.2  0.1648  0.4059    114     no     RecNo 
sigma^2.3  0.7168  0.8466    629     no    Obs_ID 

Test of Moderators (coefficients 1:6):
F(df1 = 6, df2 = 623) = 6.6538, p-val < .0001

Model Results:

              estimate      se    tval   df    pval    ci.lb   ci.ub      
Treat_modA      0.4218  0.1010  4.1766  623  <.0001   0.2235  0.6201  *** 
Treat_modAV     0.4750  0.1495  3.1778  623  0.0016   0.1815  0.7685   ** 
Treat_modAVO    0.6001  0.4574  1.3121  623  0.1900  -0.2980  1.4983      
Treat_modO      0.2655  0.2995  0.8863  623  0.3758  -0.3227  0.8537      
Treat_modOV     0.0288  0.3946  0.0729  623  0.9419  -0.7462  0.8037      
Treat_modV      0.4536  0.1108  4.0944  623  <.0001   0.2360  0.6712  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod1)*100, 2)

   R2_marginal R2_conditional 
          0.61          20.40 

Code

orchard_plot(mod1, 
             mod = "Treat_mod",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

6.2 Treatment with relevant data

Code

VCVs <- vcalc(vi = dat_short$Vd,
             cluster = dat_short$SubjectID,
             rho = 0.5)


mod1a <- rma.mv(yi = SMD,
                   V = VCVs,
                   random = list(~1|FocalSpL,
                                 ~1 | RecNo,
                                 ~1 | Obs_ID),
                   mod = ~ Treat_mod, 
                   test = "t",
                   method = "REML", 
                   sparse = TRUE,
                   data = dat_short)

# to look at A vs AV and A vs V
# TODO - make these hetero model too
summary(mod1a)

Multivariate Meta-Analysis Model (k = 589; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-834.1887  1668.3773  1680.3773  1706.6172  1680.5224   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0142  0.1190     78     no  FocalSpL 
sigma^2.2  0.1805  0.4248    102     no     RecNo 
sigma^2.3  0.7416  0.8612    589     no    Obs_ID 

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 586) = 0.0517, p-val = 0.9497

Model Results:

             estimate      se    tval   df    pval    ci.lb   ci.ub      
intrcpt        0.4249  0.1036  4.1010  586  <.0001   0.2214  0.6283  *** 
Treat_modAV    0.0542  0.1802  0.3007  586  0.7637  -0.2997  0.4081      
Treat_modV     0.0344  0.1477  0.2330  586  0.8158  -0.2557  0.3246      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

mod1b <- rma.mv(yi = SMD, 
                V = VCVs, 
                random = list(~1|FocalSpL , 
                              ~1 | RecNo, 
                              ~1 | Obs_ID), 
                mod = ~ relevel(factor(Treat_mod), ref = "AV"), 
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dat_short)

# to look at AV vs V and AV vs A (already done)
summary(mod1b)

Multivariate Meta-Analysis Model (k = 589; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-834.1887  1668.3773  1680.3773  1706.6172  1680.5224   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0142  0.1190     78     no  FocalSpL 
sigma^2.2  0.1805  0.4248    102     no     RecNo 
sigma^2.3  0.7416  0.8612    589     no    Obs_ID 

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 586) = 0.0517, p-val = 0.9497

Model Results:

                                         estimate      se     tval   df    pval 
intrcpt                                    0.4790  0.1532   3.1260  586  0.0019 
relevel(factor(Treat_mod), ref = "AV")A   -0.0542  0.1802  -0.3007  586  0.7637 
relevel(factor(Treat_mod), ref = "AV")V   -0.0198  0.1764  -0.1120  586  0.9109 
                                           ci.lb   ci.ub     
intrcpt                                   0.1781  0.7800  ** 
relevel(factor(Treat_mod), ref = "AV")A  -0.4081  0.2997     
relevel(factor(Treat_mod), ref = "AV")V  -0.3663  0.3268     

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

orchard_plot(mod1a, 
             mod = "Treat_mod",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

Code

# modeling heteroscedasticity
mod1c <- rma.mv(yi = SMD, 
                V = VCVs, 
                random = list(~1|FocalSpL , 
                              ~1 | RecNo, 
                              ~ Treat_mod | Obs_ID), 
                mod = ~ Treat_mod, 
                test = "t",
                struct = "DIAG",
                method = "REML", 
                sparse = TRUE,
                data = dat_short)

summary(mod1c)

Multivariate Meta-Analysis Model (k = 589; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-823.5686  1647.1371  1663.1371  1698.1237  1663.3867   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0126  0.1124     78     no  FocalSpL 
sigma^2.2  0.1859  0.4311    102     no     RecNo 

outer factor: Obs_ID    (nlvls = 589)
inner factor: Treat_mod (nlvls = 3)

            estim    sqrt  k.lvl  fixed  level 
tau^2.1    0.7535  0.8681    298     no      A 
tau^2.2    0.3355  0.5793     97     no     AV 
tau^2.3    0.9328  0.9658    194     no      V 

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 586) = 0.0353, p-val = 0.9653

Model Results:

             estimate      se    tval   df    pval    ci.lb   ci.ub      
intrcpt        0.4284  0.1043  4.1071  586  <.0001   0.2236  0.6333  *** 
Treat_modAV    0.0245  0.1658  0.1481  586  0.8823  -0.3010  0.3501      
Treat_modV     0.0403  0.1525  0.2642  586  0.7917  -0.2592  0.3398      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

# comparision models
anova(mod1a, mod1c)

        df       AIC       BIC      AICc    logLik     LRT   pval QE 
Full     8 1663.1371 1698.1237 1663.3867 -823.5686                NA 
Reduced  6 1680.3773 1706.6172 1680.5224 -834.1887 21.2402 <.0001 NA 

Code

orchard_plot(mod1c, 
             mod = "Treat_mod",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

6.3 Treatment with some additions

Code

# the effect of additions
# this is a part of sensitivity analysis

VCVs2 <- vcalc(vi = dat_short_add$Vd,
             cluster = dat_short_add$SubjectID,
             rho = 0.5)

mod5 <- rma.mv(yi = SMD, 
                V = VCVs2, 
                random = list(~1|FocalSpL , 
                              ~1 | RecNo, 
                              ~ Treat_mod | Obs_ID), 
                mod = ~ Treat_mod*Add_on - 1,
                test = "t",
                struct = "DIAG",
                method = "REML", 
                sparse = TRUE,
                data = dat_short_add)

summary(mod5)

Multivariate Meta-Analysis Model (k = 291; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-402.8726   805.7453   821.7453   851.0211   822.2633   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0001     38     no  FocalSpL 
sigma^2.2  0.1854  0.4305     57     no     RecNo 

outer factor: Obs_ID    (nlvls = 291)
inner factor: Treat_mod (nlvls = 2)

            estim    sqrt  k.lvl  fixed  level 
tau^2.1    0.3360  0.5796     97     no     AV 
tau^2.2    0.9062  0.9519    194     no      V 

Test of Moderators (coefficients 1:4):
F(df1 = 4, df2 = 287) = 6.0033, p-val = 0.0001

Model Results:

                      estimate      se     tval   df    pval    ci.lb   ci.ub 
Treat_modAV             0.4048  0.1417   2.8569  287  0.0046   0.1259  0.6837 
Treat_modV              0.5009  0.1283   3.9029  287  0.0001   0.2483  0.7534 
Add_onYes               0.1738  0.3505   0.4959  287  0.6204  -0.5161  0.8637 
Treat_modV:Add_onYes   -0.4974  0.4211  -1.1812  287  0.2385  -1.3264  0.3315 
                          
Treat_modAV            ** 
Treat_modV            *** 
Add_onYes                 
Treat_modV:Add_onYes      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

6.4 Treatment as an ordinal variable

Code

# testing the number of stimuli

mod4 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo, 
                             ~1 | Obs_ID), 
               mod = ~ Treat_No, 
               test = "t",
               method = "REML", 
               sparse = TRUE,
               data = dat)

summary(mod4)

Multivariate Meta-Analysis Model (k = 629; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-880.3980  1760.7961  1770.7961  1793.0008  1770.8927   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0204  0.1427     87     no  FocalSpL 
sigma^2.2  0.1520  0.3899    114     no     RecNo 
sigma^2.3  0.7135  0.8447    629     no    Obs_ID 

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 627) = 0.0263, p-val = 0.8713

Model Results:

          estimate      se    tval   df    pval    ci.lb   ci.ub    
intrcpt     0.4060  0.1722  2.3575  627  0.0187   0.0678  0.7442  * 
Treat_No    0.0209  0.1288  0.1621  627  0.8713  -0.2320  0.2737    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

bubble_plot(mod4,
             mod = "Treat_No",
             group = "RecNo",
             xlab = "The number of simuli",
             g = TRUE)

6.5 Trait type

Code

# Type of responses
mod2 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1 | FocalSpL,
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mod = ~ Type - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)

summary(mod2)

Multivariate Meta-Analysis Model (k = 629; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-872.5392  1745.0783  1757.0783  1783.7144  1757.2140   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0005     87     no  FocalSpL 
sigma^2.2  0.1378  0.3713    114     no     RecNo 
sigma^2.3  0.7145  0.8453    629     no    Obs_ID 

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 626) = 18.6766, p-val < .0001

Model Results:

                 estimate      se    tval   df    pval    ci.lb   ci.ub      
TypeBehaviour      0.5060  0.0697  7.2638  626  <.0001   0.3692  0.6428  *** 
TypeLifeHistory    0.3085  0.1195  2.5804  626  0.0101   0.0737  0.5432    * 
TypePhysiology     0.0179  0.1252  0.1426  626  0.8867  -0.2280  0.2637      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

mod2c <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL,
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mod = ~ Type,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)

summary(mod2c)

Multivariate Meta-Analysis Model (k = 629; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-872.5392  1745.0783  1757.0783  1783.7144  1757.2140   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0005     87     no  FocalSpL 
sigma^2.2  0.1378  0.3713    114     no     RecNo 
sigma^2.3  0.7145  0.8453    629     no    Obs_ID 

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 626) = 7.2108, p-val = 0.0008

Model Results:

                 estimate      se     tval   df    pval    ci.lb    ci.ub      
intrcpt            0.5060  0.0697   7.2638  626  <.0001   0.3692   0.6428  *** 
TypeLifeHistory   -0.1976  0.1240  -1.5938  626  0.1115  -0.4410   0.0459      
TypePhysiology    -0.4882  0.1304  -3.7440  626  0.0002  -0.7442  -0.2321  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

mod2d <- rma.mv(yi = SMD,
               V = VCV,
               random = list(~1|FocalSpL,
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mod = ~ relevel(Type, ref = "LifeHistory"),
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)

summary(mod2d)

Multivariate Meta-Analysis Model (k = 629; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-872.5392  1745.0783  1757.0783  1783.7144  1757.2140   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0005     87     no  FocalSpL 
sigma^2.2  0.1378  0.3713    114     no     RecNo 
sigma^2.3  0.7145  0.8453    629     no    Obs_ID 

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 626) = 7.2108, p-val = 0.0008

Model Results:

                                              estimate      se     tval   df 
intrcpt                                         0.3085  0.1195   2.5804  626 
relevel(Type, ref = "LifeHistory")Behaviour     0.1976  0.1240   1.5938  626 
relevel(Type, ref = "LifeHistory")Physiology   -0.2906  0.1546  -1.8796  626 
                                                pval    ci.lb   ci.ub    
intrcpt                                       0.0101   0.0737  0.5432  * 
relevel(Type, ref = "LifeHistory")Behaviour   0.1115  -0.0459  0.4410    
relevel(Type, ref = "LifeHistory")Physiology  0.0606  -0.5942  0.0130  . 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

orchard_plot(mod2,
             mod = "Type",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

Code

round(r2_ml(mod2)*100, 2)

   R2_marginal R2_conditional 
          3.75          19.31 

Code

# heteroscadasticity model
mod2b <- rma.mv(yi = SMD, 
               V = VCV, 
               mod = ~ Type - 1, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo, 
                             ~Type | Obs_ID), 
               struct = "DIAG",
               test = "t",
               method = "REML", 
               sparse = TRUE,
               data = dat)

summary(mod2b)

Multivariate Meta-Analysis Model (k = 629; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-853.1043  1706.2085  1722.2085  1757.7233  1722.4419   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0386  0.1964     87     no  FocalSpL 
sigma^2.2  0.1266  0.3558    114     no     RecNo 

outer factor: Obs_ID (nlvls = 629)
inner factor: Type   (nlvls = 3)

            estim    sqrt  k.lvl  fixed        level 
tau^2.1    0.6536  0.8084    412     no    Behaviour 
tau^2.2    0.3150  0.5613    112     no  LifeHistory 
tau^2.3    1.3365  1.1560    105     no   Physiology 

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 626) = 18.2414, p-val < .0001

Model Results:

                 estimate      se    tval   df    pval    ci.lb   ci.ub      
TypeBehaviour      0.5290  0.0742  7.1266  626  <.0001   0.3832  0.6748  *** 
TypeLifeHistory    0.3701  0.1065  3.4749  626  0.0005   0.1609  0.5792  *** 
TypePhysiology     0.0326  0.1525  0.2141  626  0.8305  -0.2668  0.3320      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

# make other hetero
mod2e <- rma.mv(yi = SMD,
               V = VCV,
               mod = ~ Type, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo, 
                             ~Type | Obs_ID), 
               struct = "DIAG",
               test = "t",
               method = "REML", 
               sparse = TRUE,
               data = dat)

summary(mod2e)

Multivariate Meta-Analysis Model (k = 629; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-853.1043  1706.2085  1722.2085  1757.7233  1722.4419   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0386  0.1964     87     no  FocalSpL 
sigma^2.2  0.1266  0.3558    114     no     RecNo 

outer factor: Obs_ID (nlvls = 629)
inner factor: Type   (nlvls = 3)

            estim    sqrt  k.lvl  fixed        level 
tau^2.1    0.6536  0.8084    412     no    Behaviour 
tau^2.2    0.3150  0.5613    112     no  LifeHistory 
tau^2.3    1.3365  1.1560    105     no   Physiology 

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 626) = 5.5587, p-val = 0.0040

Model Results:

                 estimate      se     tval   df    pval    ci.lb    ci.ub      
intrcpt            0.5290  0.0742   7.1266  626  <.0001   0.3832   0.6748  *** 
TypeLifeHistory   -0.1589  0.1038  -1.5317  626  0.1261  -0.3627   0.0448      
TypePhysiology    -0.4963  0.1520  -3.2646  626  0.0012  -0.7949  -0.1978   ** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

mod2f <- rma.mv(yi = SMD,
               V = VCV,
               mod = ~ relevel(Type, ref = "LifeHistory"), 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo, 
                             ~Type | Obs_ID), 
               struct = "DIAG",
               test = "t",
               method = "REML", 
               sparse = TRUE,
               data = dat)

summary(mod2f)

Multivariate Meta-Analysis Model (k = 629; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-853.1043  1706.2085  1722.2085  1757.7233  1722.4419   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0386  0.1964     87     no  FocalSpL 
sigma^2.2  0.1266  0.3558    114     no     RecNo 

outer factor: Obs_ID (nlvls = 629)
inner factor: Type   (nlvls = 3)

            estim    sqrt  k.lvl  fixed        level 
tau^2.1    0.6536  0.8084    412     no    Behaviour 
tau^2.2    0.3150  0.5613    112     no  LifeHistory 
tau^2.3    1.3365  1.1560    105     no   Physiology 

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 626) = 5.5587, p-val = 0.0040

Model Results:

                                              estimate      se     tval   df 
intrcpt                                         0.3701  0.1065   3.4749  626 
relevel(Type, ref = "LifeHistory")Behaviour     0.1589  0.1038   1.5317  626 
relevel(Type, ref = "LifeHistory")Physiology   -0.3374  0.1593  -2.1177  626 
                                                pval    ci.lb    ci.ub      
intrcpt                                       0.0005   0.1609   0.5792  *** 
relevel(Type, ref = "LifeHistory")Behaviour   0.1261  -0.0448   0.3627      
relevel(Type, ref = "LifeHistory")Physiology  0.0346  -0.6503  -0.0245    * 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

# heteroscadasticity model better than the homoscedasticity model
# note LifeHistory has lowest variation but this may be expected? 
# as it is less flexiable (e.g. the number of eggs?)
anova(mod2, mod2b)

        df       AIC       BIC      AICc    logLik     LRT   pval QE 
Full     8 1722.2085 1757.7233 1722.4419 -853.1043                NA 
Reduced  6 1757.0783 1783.7144 1757.2140 -872.5392 38.8698 <.0001 NA 

Code

orchard_plot(mod2b, 
             mod = "Type",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

6.6 trait categories

Code

# Category of responses

mod3 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo, 
                             ~1 | Obs_ID), 
               mod = ~ Category - 1, 
               test = "t",
               method = "REML", 
               sparse = TRUE,
               data = dat)

summary(mod3)

Multivariate Meta-Analysis Model (k = 629; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-866.3415  1732.6830  1754.6830  1803.4277  1755.1165   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0003     87     no  FocalSpL 
sigma^2.2  0.1411  0.3756    114     no     RecNo 
sigma^2.3  0.7176  0.8471    629     no    Obs_ID 

Test of Moderators (coefficients 1:8):
F(df1 = 8, df2 = 621) = 7.2719, p-val < .0001

Model Results:

                          estimate      se     tval   df    pval    ci.lb 
CategoryAntiPredator        0.5189  0.0964   5.3817  621  <.0001   0.3296 
CategoryCondition          -0.0169  0.1350  -0.1255  621  0.9002  -0.2820 
CategoryCostlyBehaviours    0.3963  0.1647   2.4066  621  0.0164   0.0729 
CategoryHormones            0.1743  0.2961   0.5885  621  0.5564  -0.4073 
CategoryIntake              0.6762  0.1775   3.8107  621  0.0002   0.3278 
CategoryParentalCare        0.4717  0.1113   4.2388  621  <.0001   0.2532 
CategoryPhenology           0.2005  0.1828   1.0967  621  0.2732  -0.1585 
Categoryreproduction        0.3429  0.1306   2.6251  621  0.0089   0.0864 
                           ci.ub      
CategoryAntiPredator      0.7083  *** 
CategoryCondition         0.2481      
CategoryCostlyBehaviours  0.7197    * 
CategoryHormones          0.7559      
CategoryIntake            1.0247  *** 
CategoryParentalCare      0.6903  *** 
CategoryPhenology         0.5595      
Categoryreproduction      0.5995   ** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod3)*100, 2)

   R2_marginal R2_conditional 
          4.28          20.01 

Code

orchard_plot(mod3, 
             mod = "Category",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             angle = 45,
             branch.size = 3)

6.7 Predactor guild

Code

# Predactor guild
# quite heterogeneous
# TODO this could be in random effects - think abou thtis a bit later
mod6 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ PredGuild - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)

summary(mod6)

Multivariate Meta-Analysis Model (k = 629; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-873.0918  1746.1836  1764.1836  1804.0945  1764.4772   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0272  0.1648     87     no  FocalSpL 
sigma^2.2  0.1504  0.3878    114     no     RecNo 
sigma^2.3  0.7153  0.8458    629     no    Obs_ID 

Test of Moderators (coefficients 1:6):
F(df1 = 6, df2 = 623) = 7.2101, p-val < .0001

Model Results:

                  estimate      se     tval   df    pval    ci.lb   ci.ub      
PredGuildBird       0.4945  0.0815   6.0658  623  <.0001   0.3344  0.6546  *** 
PredGuildFish      -0.1862  0.9303  -0.2002  623  0.8414  -2.0131  1.6407      
PredGuildMammal     0.3053  0.1335   2.2874  623  0.0225   0.0432  0.5674    * 
PredGuildmixed      0.1110  0.2573   0.4313  623  0.6664  -0.3943  0.6162      
PredGuildNS         0.6943  0.3095   2.2437  623  0.0252   0.0866  1.3021    * 
PredGuildReptile    0.4022  0.3691   1.0896  623  0.2763  -0.3226  1.1270      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod6)*100, 2)

   R2_marginal R2_conditional 
          2.22          21.67 

Code

orchard_plot(mod6, 
             mod = "PredGuild",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

6.8 Setting

Code

# Setting

mod7 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ Setting - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)

summary(mod7)

Multivariate Meta-Analysis Model (k = 629; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-878.0016  1756.0032  1768.0032  1794.6393  1768.1389   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0239  0.1547     87     no  FocalSpL 
sigma^2.2  0.1445  0.3801    114     no     RecNo 
sigma^2.3  0.7129  0.8444    629     no    Obs_ID 

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 626) = 13.8338, p-val < .0001

Model Results:

                     estimate      se     tval   df    pval    ci.lb   ci.ub 
SettingField           0.4528  0.0764   5.9271  626  <.0001   0.3028  0.6028 
SettingLab             0.4251  0.1561   2.7230  626  0.0067   0.1185  0.7316 
SettingSemi-natural   -0.1997  0.4319  -0.4622  626  0.6441  -1.0479  0.6486 
                         
SettingField         *** 
SettingLab            ** 
SettingSemi-natural      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod7)*100, 2)

   R2_marginal R2_conditional 
          0.68          19.66 

Code

orchard_plot(mod7, 
             mod = "Setting",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

6.9 Season

Code

# Season

mod8 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ Season - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)

mod8b <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ Season,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)

summary(mod8b)

Multivariate Meta-Analysis Model (k = 629; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-879.5540  1759.1080  1769.1080  1791.3128  1769.2047   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0287  0.1695     87     no  FocalSpL 
sigma^2.2  0.1444  0.3800    114     no     RecNo 
sigma^2.3  0.7129  0.8444    629     no    Obs_ID 

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 627) = 1.4957, p-val = 0.2218

Model Results:

                    estimate      se    tval   df    pval    ci.lb   ci.ub      
intrcpt               0.3858  0.0823  4.6887  627  <.0001   0.2242  0.5473  *** 
Seasonnon-breeding    0.1686  0.1379  1.2230  627  0.2218  -0.1021  0.4394      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod8)*100, 2)

   R2_marginal R2_conditional 
          0.65          20.06 

Code

orchard_plot(mod8,
             mod = "Season",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

6.10 Design

Code

# Design
mod9 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ Design - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)
summary(mod9)

Multivariate Meta-Analysis Model (k = 629; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-880.3002  1760.6003  1770.6003  1792.8050  1770.6969   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0180  0.1343     87     no  FocalSpL 
sigma^2.2  0.1536  0.3919    114     no     RecNo 
sigma^2.3  0.7140  0.8450    629     no    Obs_ID 

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 627) = 19.5259, p-val < .0001

Model Results:

              estimate      se    tval   df    pval   ci.lb   ci.ub      
Designamong     0.4057  0.0893  4.5458  627  <.0001  0.2305  0.5810  *** 
Designwithin    0.4538  0.0899  5.0472  627  <.0001  0.2772  0.6303  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

mod9b <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ Design,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)

summary(mod9b)

Multivariate Meta-Analysis Model (k = 629; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-880.3002  1760.6003  1770.6003  1792.8050  1770.6969   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0180  0.1343     87     no  FocalSpL 
sigma^2.2  0.1536  0.3919    114     no     RecNo 
sigma^2.3  0.7140  0.8450    629     no    Obs_ID 

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 627) = 0.1760, p-val = 0.6750

Model Results:

              estimate      se    tval   df    pval    ci.lb   ci.ub      
intrcpt         0.4057  0.0893  4.5458  627  <.0001   0.2305  0.5810  *** 
Designwithin    0.0480  0.1145  0.4195  627  0.6750  -0.1768  0.2729      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod9)*100, 2)

   R2_marginal R2_conditional 
          0.06          19.43 

Code

orchard_plot(mod9,
             mod = "Design",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

6.11 Response period

Code

# Response period
mod10 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ ResponsePeriod - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)
summary(mod10)

Multivariate Meta-Analysis Model (k = 629; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-879.6272  1759.2544  1769.2544  1791.4591  1769.3510   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0207  0.1439     87     no  FocalSpL 
sigma^2.2  0.1496  0.3868    114     no     RecNo 
sigma^2.3  0.7142  0.8451    629     no    Obs_ID 

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 627) = 20.2188, p-val < .0001

Model Results:

                      estimate      se    tval   df    pval   ci.lb   ci.ub 
ResponsePeriodafter     0.3686  0.0873  4.2215  627  <.0001  0.1972  0.5401 
ResponsePeriodduring    0.5106  0.0959  5.3257  627  <.0001  0.3223  0.6989 
                          
ResponsePeriodafter   *** 
ResponsePeriodduring  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod10)*100, 2)

   R2_marginal R2_conditional 
          0.56          19.71 

Code

orchard_plot(mod10,
             mod = "ResponsePeriod",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

6.12 Control type

Code

# control type
mod11 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ ControlType - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)
summary(mod11)

Multivariate Meta-Analysis Model (k = 629; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-878.4467  1756.8935  1770.8935  1801.9577  1771.0750   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0191  0.1381     87     no  FocalSpL 
sigma^2.2  0.1539  0.3922    114     no     RecNo 
sigma^2.3  0.7161  0.8462    629     no    Obs_ID 

Test of Moderators (coefficients 1:4):
F(df1 = 4, df2 = 625) = 9.7230, p-val < .0001

Model Results:

                        estimate      se    tval   df    pval    ci.lb   ci.ub 
ControlTypeBlank          0.4485  0.0970  4.6246  625  <.0001   0.2580  0.6389 
ControlTypedisturbance    0.4477  0.1479  3.0274  625  0.0026   0.1573  0.7381 
ControlTypemix            0.4547  0.9640  0.4717  625  0.6373  -1.4384  2.3478 
ControlTypeNonPred        0.4147  0.0853  4.8642  625  <.0001   0.2473  0.5821 
                            
ControlTypeBlank        *** 
ControlTypedisturbance   ** 
ControlTypemix              
ControlTypeNonPred      *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod11)*100, 2)

   R2_marginal R2_conditional 
          0.03          19.48 

Code

orchard_plot(mod11,
             mod = "ControlType",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

6.13 Sex

Code

# sex
# TODO - this could be interesting
# what is in males and females
mod12 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ Sex - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)
summary(mod12)

Multivariate Meta-Analysis Model (k = 629; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-876.9390  1753.8780  1765.8780  1792.5141  1766.0137   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0249  0.1577     87     no  FocalSpL 
sigma^2.2  0.1323  0.3637    114     no     RecNo 
sigma^2.3  0.7147  0.8454    629     no    Obs_ID 

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 626) = 15.0953, p-val < .0001

Model Results:

         estimate      se    tval   df    pval   ci.lb   ci.ub      
Sexboth    0.5117  0.0806  6.3506  626  <.0001  0.3535  0.6700  *** 
SexF       0.2536  0.1061  2.3904  626  0.0171  0.0453  0.4620    * 
SexM       0.4141  0.1356  3.0546  626  0.0023  0.1479  0.6804   ** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod12)*100, 2)

   R2_marginal R2_conditional 
          1.32          19.11 

Code

orchard_plot(mod12,
             mod = "Sex",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

Code

# shoter data for just males and females
# hetero but no sex effect here
dat_sex <- dat %>% filter(Sex != "both")

VCV3 <- vcalc(vi = dat_sex$Vd,
             cluster = dat_sex$SubjectID,
             rho = 0.5)

mod12a <- rma.mv(yi = SMD, 
                 V = VCV3, 
                 mod = ~ Sex, 
                 random = list(~1|FocalSpL , 
                               ~1 | RecNo, 
                               ~1 | Obs_ID), 
                 #struct = "DIAG",
                 test = "t",
                 method = "REML", 
                 sparse = TRUE,
                 data = dat_sex)

mod12b <- rma.mv(yi = SMD, 
                V = VCV3, 
                mod = ~ Sex, 
                random = list(~1|FocalSpL , 
                              ~1 | RecNo, 
                              ~Sex | Obs_ID), 
                struct = "DIAG",
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dat_sex)

summary(mod12b)

Multivariate Meta-Analysis Model (k = 256; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-285.0946   570.1892   582.1892   603.4132   582.5293   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0931  0.3052     35     no  FocalSpL 
sigma^2.2  0.0000  0.0009     48     no     RecNo 

outer factor: Obs_ID (nlvls = 256)
inner factor: Sex    (nlvls = 2)

            estim    sqrt  k.lvl  fixed  level 
tau^2.1    0.3017  0.5493    155     no      F 
tau^2.2    0.7744  0.8800    101     no      M 

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 254) = 2.5399, p-val = 0.1122

Model Results:

         estimate      se     tval   df    pval    ci.lb   ci.ub      
intrcpt    0.4627  0.1327   3.4863  254  0.0006   0.2013  0.7240  *** 
SexF      -0.2068  0.1298  -1.5937  254  0.1122  -0.4624  0.0487      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

anova(mod12a, mod12b)

        df      AIC      BIC     AICc    logLik     LRT   pval QE 
Full     6 582.1892 603.4132 582.5293 -285.0946                NA 
Reduced  5 598.6866 616.3733 598.9286 -294.3433 18.4974 <.0001 NA 

Code

orchard_plot(mod12b,
             mod = "Sex",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

6.14 Age

Code

# age
mod13 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ Age - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)
summary(mod13)

Multivariate Meta-Analysis Model (k = 629; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-877.7166  1755.4331  1769.4331  1800.4974  1769.6146   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0191  0.1383     87     no  FocalSpL 
sigma^2.2  0.1534  0.3917    114     no     RecNo 
sigma^2.3  0.7165  0.8464    629     no    Obs_ID 

Test of Moderators (coefficients 1:4):
F(df1 = 4, df2 = 625) = 9.8896, p-val < .0001

Model Results:

      estimate      se    tval   df    pval    ci.lb   ci.ub      
AgeA    0.4465  0.0717  6.2273  625  <.0001   0.3057  0.5873  *** 
AgeE    0.3414  0.3552  0.9611  625  0.3369  -0.3562  1.0389      
AgeJ    0.3069  0.3553  0.8639  625  0.3880  -0.3908  1.0047      
AgeN    0.2993  0.1998  1.4979  625  0.1347  -0.0931  0.6917      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod13)*100, 2)

   R2_marginal R2_conditional 
          0.30          19.65 

Code

orchard_plot(mod13,
             mod = "Age",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

6.15 Predactor type

Code

# type of predator

dat$PredTo[dat$PredTo == ""] <- NA
mod14 <- rma.mv(yi = SMD, 
               V = VCV, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ PredTo - 1,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)
summary(mod14)

Multivariate Meta-Analysis Model (k = 610; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-848.5688  1697.1375  1709.1375  1735.5887  1709.2775   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0038  0.0615     82     no  FocalSpL 
sigma^2.2  0.1663  0.4078    108     no     RecNo 
sigma^2.3  0.7082  0.8416    610     no    Obs_ID 

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 607) = 11.8741, p-val < .0001

Model Results:

         estimate      se    tval   df    pval    ci.lb   ci.ub      
PredToA    0.4386  0.0848  5.1721  607  <.0001   0.2720  0.6051  *** 
PredToB    0.3171  0.3581  0.8856  607  0.3762  -0.3861  1.0203      
PredToN    0.3581  0.0960  3.7308  607  0.0002   0.1696  0.5466  *** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod14)*100, 2)

   R2_marginal R2_conditional 
          0.20          19.53 

Code

orchard_plot(mod14,
             mod = "PredTo",
             group = "RecNo",
             xlab = "Standardised mean differnece (SMD)")

6.16 Treatment duration

Code

# treatment duration

dat$ln_duration <- log(dat$duration_days)

mod15 <- rma.mv(yi = SMD,
               V = VCV,
               random = list(~1|FocalSpL,
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ ln_duration,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)
summary(mod15)

Multivariate Meta-Analysis Model (k = 629; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-875.1616  1750.3231  1760.3231  1782.5278  1760.4197   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0001     87     no  FocalSpL 
sigma^2.2  0.1592  0.3989    114     no     RecNo 
sigma^2.3  0.7139  0.8449    629     no    Obs_ID 

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 627) = 10.5064, p-val = 0.0013

Model Results:

             estimate      se     tval   df    pval    ci.lb    ci.ub      
intrcpt        0.3799  0.0652   5.8293  627  <.0001   0.2519   0.5078  *** 
ln_duration   -0.0515  0.0159  -3.2414  627  0.0013  -0.0827  -0.0203   ** 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

mod16 <- rma.mv(yi = SMD,
               V = VCV,
               random = list(~1|FocalSpL,
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ ln_duration*Type,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat)
summary(mod16)

Multivariate Meta-Analysis Model (k = 629; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-864.7519  1729.5039  1747.5039  1787.4148  1747.7975   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0001     87     no  FocalSpL 
sigma^2.2  0.1491  0.3861    114     no     RecNo 
sigma^2.3  0.7086  0.8418    629     no    Obs_ID 

Test of Moderators (coefficients 2:6):
F(df1 = 5, df2 = 623) = 4.5323, p-val = 0.0005

Model Results:

                             estimate      se     tval   df    pval    ci.lb 
intrcpt                        0.4564  0.0732   6.2347  623  <.0001   0.3126 
ln_duration                   -0.0481  0.0174  -2.7676  623  0.0058  -0.0822 
TypeLifeHistory               -0.3590  0.3642  -0.9858  623  0.3246  -1.0742 
TypePhysiology                -0.5372  0.1648  -3.2604  623  0.0012  -0.8607 
ln_duration:TypeLifeHistory    0.0913  0.1027   0.8893  623  0.3742  -0.1103 
ln_duration:TypePhysiology     0.0664  0.0514   1.2925  623  0.1967  -0.0345 
                               ci.ub      
intrcpt                       0.6001  *** 
ln_duration                  -0.0140   ** 
TypeLifeHistory               0.3562      
TypePhysiology               -0.2136   ** 
ln_duration:TypeLifeHistory   0.2929      
ln_duration:TypePhysiology    0.1673      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod15)*100, 2)

   R2_marginal R2_conditional 
          4.25          21.70 

Code

bubble_plot(mod15,
             mod = "ln_duration",
             group = "RecNo",
             xlab = "log(duration in days)",
             g = TRUE) +
    geom_point(data = dat,
    aes(x = ln_duration, y = SMD,
    color = Type,
    fill = Type,
    size = 1/sqrt(Vd)), alpha = 0.6) +
    scale_color_discrete() + #+ # how to put the legend for colour
    guides(color = "legend")

Code

#p + geom_point(aes(colour = Type))
#scale_colour_manual(values = c("red", "blue", "green"))

7 Meta-regression: multi-moderator (mostly)

7.1 full models

Code

#######################
# Mulit-variable models
#######################

dat_short$sln_duration <- scale(log(dat_short$duration_days))

mod_full <- rma.mv(yi = SMD, 
               V = VCVs, 
               random = list(~1|FocalSpL , 
                             ~1 | RecNo,
                             ~1 | Obs_ID),
               mods =  ~ #Design +
                         sln_duration*Type +
                         sln_duration*Treat_mod +
                         Sex,
               test = "t",
               method = "REML",
               sparse = TRUE,
               data = dat_short)
summary(mod_full)

Multivariate Meta-Analysis Model (k = 589; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-809.4704  1618.9407  1648.9407  1714.3083  1649.7963   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0003     78     no  FocalSpL 
sigma^2.2  0.1580  0.3975    102     no     RecNo 
sigma^2.3  0.7361  0.8579    589     no    Obs_ID 

Test of Moderators (coefficients 2:12):
F(df1 = 11, df2 = 577) = 2.6199, p-val = 0.0029

Model Results:

                              estimate      se     tval   df    pval    ci.lb 
intrcpt                         0.5125  0.1106   4.6348  577  <.0001   0.2953 
sln_duration                   -0.1431  0.0912  -1.5694  577  0.1171  -0.3223 
TypeLifeHistory                -0.0395  0.4523  -0.0872  577  0.9305  -0.9279 
TypePhysiology                 -0.5792  0.1791  -3.2340  577  0.0013  -0.9309 
Treat_modAV                     0.1811  0.1795   1.0093  577  0.3132  -0.1713 
Treat_modV                      0.0344  0.1511   0.2276  577  0.8201  -0.2624 
SexF                           -0.3166  0.1340  -2.3629  577  0.0185  -0.5799 
SexM                           -0.1304  0.1570  -0.8304  577  0.4066  -0.4387 
sln_duration:TypeLifeHistory    0.1825  0.4565   0.3997  577  0.6895  -0.7142 
sln_duration:TypePhysiology     0.3173  0.2080   1.5252  577  0.1278  -0.0913 
sln_duration:Treat_modAV       -0.2158  0.1955  -1.1039  577  0.2701  -0.5997 
sln_duration:Treat_modV        -0.0722  0.1436  -0.5028  577  0.6153  -0.3541 
                                ci.ub      
intrcpt                        0.7297  *** 
sln_duration                   0.0360      
TypeLifeHistory                0.8490      
TypePhysiology                -0.2274   ** 
Treat_modAV                    0.5336      
Treat_modV                     0.3311      
SexF                          -0.0534    * 
SexM                           0.1780      
sln_duration:TypeLifeHistory   1.0792      
sln_duration:TypePhysiology    0.7258      
sln_duration:Treat_modAV       0.1682      
sln_duration:Treat_modV        0.2098      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

round(r2_ml(mod_full)*100, 2)

   R2_marginal R2_conditional 
          9.04          25.12 

Code

orchard_plot(mod_full,
             mod = "Type",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

Code

orchard_plot(mod_full,
             mod = "Treat_mod",
             group = "RecNo", 
             xlab = "Standardised mean differnece (SMD)",
             branch.size = 3)

Code

int_type <- mod_results(mod_full, mod = "sln_duration", group = "RecNo", weights = "prop",
                                   by = "Type")

bubble_plot(int_type, group = "RecNo", mod = "sln_duration", xlab = "ln(duration in days)",
                     legend.pos = "top.left", condition.nrow = 3)

Code

int_trt <- mod_results(mod_full, mod = "sln_duration", group = "RecNo", weights = "prop",
                        by = "Treat_mod")

bubble_plot(int_trt, group = "RecNo", mod = "sln_duration", xlab = "ln(duration in days)",
            legend.pos = "top.left", condition.nrow = 3)

Code

# mulit-model selection
candidates <- dredge(mod_full, trace = 2)

  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=                                                                     |   2%
  |                                                                            
  |==                                                                    |   3%
  |                                                                            
  |===                                                                   |   5%
  |                                                                            
  |====                                                                  |   6%
  |                                                                            
  |=====                                                                 |   8%
  |                                                                            
  |=======                                                               |   9%
  |                                                                            
  |========                                                              |  11%
  |                                                                            
  |=========                                                             |  12%
  |                                                                            
  |==========                                                            |  14%
  |                                                                            
  |===========                                                           |  16%
  |                                                                            
  |============                                                          |  17%
  |                                                                            
  |=============                                                         |  19%
  |                                                                            
  |==============                                                        |  20%
  |                                                                            
  |===============                                                       |  22%
  |                                                                            
  |================                                                      |  23%
  |                                                                            
  |========================                                              |  34%
  |                                                                            
  |=========================                                             |  36%
  |                                                                            
  |=================================                                     |  47%
  |                                                                            
  |==================================                                    |  48%
  |                                                                            
  |==============================================                        |  66%
  |                                                                            
  |===============================================                       |  67%
  |                                                                            
  |==================================================                    |  72%
  |                                                                            
  |===================================================                   |  73%
  |                                                                            
  |====================================================================  |  97%
  |                                                                            
  |===================================================================== |  98%

Code

# displays delta AICc <2
candidates_aic2 <- subset(candidates, delta < 5) 
# model averaging
mr_averaged_aic2 <- summary(model.avg(candidates, delta < 5)) 

# relative importance of each predictor for all the models
importance <- sw(candidates)

8 Publication bias

8.1 Funnel plot: uni-moderator

8.2 Funnel plot: multi-moderator

Code

funnel(mod0r, 
       yaxis="seinv",
       type = "rstudent")

Code

funnel(mod_full, 
       yaxis="seinv",
       type = "rstudent")

8.3 Egger regression: uni-moderator

Code

# Egger

dat$effectN <- (dat$Ncontrol * dat$NTreat) / (dat$Ncontrol + dat$NTreat)
dat$sqeffectN <- sqrt(dat$effectN)

mod0e <- rma.mv(yi = SMD,
                V = VCV,
                mods = ~ sqeffectN,
                random = list(#~1 | Phylo,
                  ~1 | FocalSpL,
                  ~1 | RecNo,
                  #~1 | SubjectID, # incoprated as VCV
                  ~1 | Obs_ID),
                #R = list(Phylo = cor_tree),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dat)

summary(mod0e)

Multivariate Meta-Analysis Model (k = 629; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-880.0117  1760.0234  1770.0234  1792.2282  1770.1200   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0112  0.1059     87     no  FocalSpL 
sigma^2.2  0.1512  0.3888    114     no     RecNo 
sigma^2.3  0.7182  0.8475    629     no    Obs_ID 

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 627) = 0.9653, p-val = 0.3262

Model Results:

           estimate      se     tval   df    pval    ci.lb   ci.ub      
intrcpt      0.5260  0.1248   4.2163  627  <.0001   0.2810  0.7710  *** 
sqeffectN   -0.0287  0.0292  -0.9825  627  0.3262  -0.0861  0.0287      

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

bubble_plot(mod0e,
            mod = "sqeffectN",
            group = "RecNo",
            xlab = "Effective N",
            g = TRUE)

8.4 Decline effect: uni-moderator

Code

# decline effect
mod0d <- rma.mv(yi = SMD,
                V = VCV,
                mods = ~ Year,
                random = list(#~1 | Phylo,
                  ~1 | FocalSpL,
                  ~1 | RecNo,
                  #~1 | SubjectID, # incoprated as VCV
                  ~1 | Obs_ID),
                #R = list(Phylo = cor_tree),
                test = "t",
                method = "REML", 
                sparse = TRUE,
                data = dat)

summary(mod0d)

Multivariate Meta-Analysis Model (k = 629; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-880.0555  1760.1110  1770.1110  1792.3157  1770.2076   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0205  0.1431     87     no  FocalSpL 
sigma^2.2  0.1515  0.3892    114     no     RecNo 
sigma^2.3  0.7145  0.8453    629     no    Obs_ID 

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 627) = 0.4339, p-val = 0.5103

Model Results:

         estimate       se     tval   df    pval     ci.lb    ci.ub    
intrcpt   12.8756  18.8905   0.6816  627  0.4957  -24.2206  49.9719    
Year      -0.0062   0.0094  -0.6587  627  0.5103   -0.0246   0.0123    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

bubble_plot(mod0d,
            mod = "Year",
            group = "RecNo",
            xlab = "Publication year",
            g = TRUE)

8.5 All together

Code

# full model
dat_short$effectN <- (dat_short$Ncontrol * dat_short$NTreat) / (dat_short$Ncontrol + dat_short$NTreat)
dat_short$sqeffectN <- sqrt(dat_short$effectN)

mod_fulle <- rma.mv(yi = SMD, 
                   V = VCVs, 
                   random = list(~1|FocalSpL , 
                                 ~1 | RecNo,
                                 ~1 | Obs_ID),
                   mods =  ~ sqeffectN +
                     Year +
                     sln_duration*Type +
                     sln_duration*Treat_mod +
                     Sex,
                   test = "t",
                   method = "REML",
                   sparse = TRUE,
                   data = dat_short)
summary(mod_fulle)

Multivariate Meta-Analysis Model (k = 589; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-806.8644  1613.7289  1647.7289  1721.7532  1648.8276   

Variance Components:

            estim    sqrt  nlvls  fixed    factor 
sigma^2.1  0.0000  0.0003     78     no  FocalSpL 
sigma^2.2  0.1615  0.4019    102     no     RecNo 
sigma^2.3  0.7383  0.8592    589     no    Obs_ID 

Test of Moderators (coefficients 2:14):
F(df1 = 13, df2 = 575) = 2.2130, p-val = 0.0081

Model Results:

                              estimate       se     tval   df    pval     ci.lb 
intrcpt                         0.2443  22.9599   0.0106  575  0.9915  -44.8511 
sqeffectN                      -0.0109   0.0338  -0.3213  575  0.7481   -0.0773 
Year                            0.0002   0.0114   0.0135  575  0.9892   -0.0223 
sln_duration                   -0.1351   0.0951  -1.4200  575  0.1561   -0.3219 
TypeLifeHistory                -0.0249   0.4552  -0.0548  575  0.9563   -0.9189 
TypePhysiology                 -0.5791   0.1846  -3.1369  575  0.0018   -0.9417 
Treat_modAV                     0.1784   0.1838   0.9708  575  0.3321   -0.1826 
Treat_modV                      0.0332   0.1534   0.2163  575  0.8288   -0.2682 
SexF                           -0.3173   0.1345  -2.3586  575  0.0187   -0.5815 
SexM                           -0.1342   0.1587  -0.8456  575  0.3982   -0.4458 
sln_duration:TypeLifeHistory    0.1662   0.4601   0.3613  575  0.7180   -0.7374 
sln_duration:TypePhysiology     0.3149   0.2150   1.4641  575  0.1437   -0.1075 
sln_duration:Treat_modAV       -0.2187   0.1997  -1.0949  575  0.2740   -0.6110 
sln_duration:Treat_modV        -0.0787   0.1453  -0.5415  575  0.5884   -0.3642 
                                ci.ub     
intrcpt                       45.3398     
sqeffectN                      0.0556     
Year                           0.0226     
sln_duration                   0.0518     
TypeLifeHistory                0.8691     
TypePhysiology                -0.2165  ** 
Treat_modAV                    0.5395     
Treat_modV                     0.3346     
SexF                          -0.0531   * 
SexM                           0.1775     
sln_duration:TypeLifeHistory   1.0698     
sln_duration:TypePhysiology    0.7372     
sln_duration:Treat_modAV       0.1736     
sln_duration:Treat_modV        0.2068     

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Code

dat_fulle <- qdrg(object = mod_fulle, 
                  data = dat_short)
# marginalized overall mean at vi = 0 and year.c = 0; also weights = "prop" or "cells" average things over proportionally. if not specified, all groups (levels) get the same weights
# res_fulle1 <- emmeans(dat_fulle, 
#                      specs = ~ sqeffectN,
#                      df = mod_fulle$ddf, 
#                      weights = "prop")